Equation of state of hard oblate ellipsoids by replica exchange Monte Carlo 
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We implemented the replica exchange Monte Carlo technique to produce the equation of state of 
hard 1:5 aspect-ratio oblate ellipsoids for a wide density range. For this purpose, we considered the 
analytical approximation of the overlap distance given by Bern and Pechukas and the exact numerical 
solution given by Perram and Wertheim. For both cases we capture the expected isotropic-nematic 
transition at low densities and a nematic-crystal transition at larger densities. For the exact case, 
these transitions occur at the volume fraction 0.341, and in the interval 0.584 — 0.605, respectively. 



Anisotropic molecules and colloids show a strong ten- 
dency to assemble into complex structures which range 
from liquid-crystals [l|, Q to empty liquids [sl, 0] • There 
are basically two approaches to model these systems: site 
by site or atomistic, and coarse graining [l|,[2|. In the first 
approach detail is usually gained at the expenses of a 
larger computational cost. In the second, many sites are 
grouped into entities which, in many cases, do not show a 
spherical symmetry. In such cases, ellipsoidal pair poten- 
tials have been widely used [El-Q • Among them, probably 
the most popular is the soft Gay-Berne interaction |8Ml2] , 
which is based on the hard core overlap distance intro- 
duced by Berne and Pechukas (BP) [l^. 

An ellipsoid can be handled to match the shape of 
many molecules and colloids. For instance, the hard core 
of laponites, a clay with a well-defined disk-like shape, 
can be approached by oblates having a 1:25 aspect- 
ratio 0, 4, 14-16]. Then, charge or Yukawa sites can 
be added to match their electrostatic properties. This 
surely leads to a considerable decrease of computational 
cost since otherwise the clay body must be modeled with 
a large number of beads to mimic its hard core [14]. For 
this purpose, however, it is desirable to know the exact 
shape of the model hard core, as well as its corresponding 
volume and surface. Unfortunately, the BP hard poten- 
tial does not provide a defined shape and volume. This 
may explain the subsequent efforts to efficiently solve the 
exact ellipsoidal hard core interaction (El, Q . In this re- 
gard, Paramonov and Yalirakia W , based on the work of 
Perram and Wertheim (PW) [1, Q , recently presented a 
numerical algorithm for determining the directional con- 
tact distance of two generic ellipsoids, which can be used 
to detect overlaps. 

The BP hard potential |l^ is analytical, mathemat- 
ically simple, easy to implement, fast to compute, and 
it can be used to study the condensed phase of pro- 
late or oblate particles via simulations (ol-ITTl IT7| . In 
this approach molecules are represented with an uni- 
axial ellipsoidal Gaussian and their interaction is then 
related to their overlap integral. This way, a coarse- 
grained potential is built which successfully captures 
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the anisotropic nature of the entities. The expression 
for the distance between the geometric centers of the 
ellipsoids when the particles are at contact, cr^p, is 
given by asp = cr±(l - IxiA'^^^ + A*^"^])"^/^, where 
A'^^^ = {vji • Ui ± Tji • Uj)^/(1 ± X^i • Uj), Ui and u_^- are 
the versors along the axial axis of each particle, and r is 
the versor along the line joining the geometric centers. 
Here, x = {^f\ ~ ^i.)/{^f\ + is the anisotropy param- 
eter, where a\\ and a_\_ are the parallel and perpendicular 
diameters with respect to the axial axis, respectively. 

On the other hand, an exact solution for determining 
whether two ellipsoids overlap or not can be numerically 
obtained. An ellipsoidal surface centered at and ori- 
ented according to is given by the following quadratic 
form Ai{re) = (rg — ) ^ • • (rg — ) , where is a point 
at the surface A(re) = 1, = U^(ui) • • U(ui), 
U(ui) is the rotation matrix, is its transpose, and 
D = ^ ^ . aiii (g) Ci {e-i is a versor along the particle prin- 
cipal axis). Let's consider two arbitrarily oriented ellip- 
soids i and j at contact at point Vc. The vector normal to 
the i surface at Vc is n^(rc) = A^ • (fc — r^). A similar rela- 
tion can be written for the vector normal to the j surface 
Uj. Since the tangent plane is common for both ellip- 
soids, the normal versors hi = ni/\ni\ and hj = nj/|nj| 
fulffil hi^Tc) + nj(rc) = 0- This condition was origi- 
nally employed by PW for developing an algorithm to 
numerically determine the point Vc [H, Q . In their work 
the Elliptic Contact Function (EOF) is introduced. This 
function contains the above given information and al- 
lows to determine the distance of the closest approach. 
Later, the ECF procedure was reviewed by Paramonov 
and Yaliraki, who contributed with a clear geometrical 
interpretation of the PW approach [7]. In particular, 
the expression for the function that connects the parti- 
cles centers trough the geometrical place where the vec- 
tors V^i(xc) and V^j(xc) are antiparallel is given by 0] 
Xc(A) = {XAi + (1 - X)Aj)-^ ■ {XAi • + (1 - A)A^- • r^-), 
where AG [0, 1] is a scalar parameter. Note that for A = 1 
and the geometric centers of the ellipsoids i and j are 
obtained, respectively. The contact point Vc lies on this 
trajectory and corresponds to a unique value of A, Ac, 
which fuffills Ac G (0,1). Furthermore, Ai{rc) = Aj{rc), 
with Vc = Xc(Ac). 

With the above expressions it is easy to implement 
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an iterative procedure to yield Vc with the desired pre- 
cision. In particular, we start by evaluating A (A) = 
yli(xc(A)) - ^^•(xc(A)) for A = 0.5. A positive A(A) 
means 1 > Ac > A and so, we increase A in such a way to 
reduce in half the interval. Conversely, a negative A (A) 
means < Ac < A and we accordingly decrease A. This 
way the interval is reduced as 1/2^, being n the number 
of iterations. Approximately 20 iterations yield an error 
of A smaller that 1 x 10~^. Note that the involved oper- 
ations are products and summations which translate into 
a relatively fast computation. 

Additionally, the contact parameter Ac is the extreme 
value of the linear combination of the quadratic forms 
A(xc) and ^j(xc), i. e., S{X) = AA(xc(A)) + (1 - 
A)^^(xc(A)) and < 5(A) < 5(Ac) [3]. This property 
can also be used to numerically determine Ac- The con- 
tact parameter defines the PW contact distance, crpw = 
r/ y^S{Xc) [5|-01- Consequently, ellipsoids having their 
geometric centers separated at a distance r < apw over- 
lap whereas they do not for r > apw In particular, the 
BP analytical expression for the contact distance corre- 
sponds to a BP = r/ ^/S{l/2) [7]. This expression makes 
clear that cbp ^ ^pw- 

Even though the analytical expressions for determining 
whether or not two ellipsoids overlap are relatively fast 
to compute, sampling from crowded systems is always a 
difficult task [18]. Thus, we implemented the replica ex- 
change Monte Carlo methodology, which is well proven 
to assist the systems to reach equilibrium at difficult 
(high density / low temperature) conditions [9, 19-23]. 
Since we are dealing with hard ellipsoids we must perform 
the replica expansion in pressure. Hence, the partition 
function of the extended ensemble is given by [2^ 
Qextended = YY^LiQNTPi, whcrc QNTPi is the partition 
function of the isobaric-isothermal ensemble of the sys- 
tem at pressure P^, temperature T, and particle number 
N . This extended ensemble is sampled by combining 
standard NTPi simulations on each replica and swap 
moves at the replica level. These swap moves are per- 
formed by means of the following acceptance probabil- 
ity [23] i^cc = min(l,exp[/3(i^ - Pj){V, - Vj)]), where 
Vi — Vj is the volume difference between replicas i and j. 
More details on this method are given in refs. [HI, [23^ . 

Simulations are started by randomly placing and ori- 
enting the ellipsoids (avoiding overlaps), so that the ini- 
tial volume fraction is cp = Vep = 0.2, where p is the 
number density, = 47r(j||<j^/3 is the ellipsoid volume 
(for both studied models), a± = Scry, and cry is taken as 
the length unit. We first perform 2 x 10^^ trial moves 
at the desired state points, during which we observe the 
replicas reaching a stationary state. We then sample by 
performing additional 2 x 10^^ trials. This work is per- 
formed by considering = 100 and = 64. 

The exact hard 1:5 oblate ellipsoidal model can be 
studied since, on the one hand, the exact iterative proce- 
dure to solve the apw overlap distance is relatively fast, 
and on the other, the numerical solution is computed 
only for those cases where r < app (otherwise ellipsoids 
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FIG. 1. a) Probability distribution functions (PDFs) of vol- 
ume fraction fluctuations for each of the rir = 64 pressure 
values, and for the exact hard 1:5 oblate ellipsoidal model, 
b) Equations of state, Z{(p). c) Isothermal compressibilities, 
x(^), obtained from density fluctuations, d) Order parame- 
ters, Q6{(p). For panels b), c), and d), dark circles and light 
squares correspond to the exact PW and the BP analytical 
solution for the overlap distance, respectively. Vertical dotted 
lines highlight the PW transitions. 



do not overlap). The probability distribution functions 
(PDFs) obtained for this system are shown in Fig. [T]a). 
There are 64 curves corresponding to each fixed pressure, 
which increases from left to right. The general trend of 
the PDFs is to get narrower and higher with increasing 
pressure evidencing a decrease of the isothermal com- 
pressibility X = Sp/S{/3P). However, dit (p ^ 0.34 and at 
(p ^ 0.60 this trend is disrupted. Here the PDFs turn 
wider, shorter, and distorted, pointing out phase tran- 
sitions. In particular, at (p ^ 0.60 PDFs are bimodal, 
which is the typical behavior of a first order transition (a 
discontinuity of the pressure as a function of the density). 
Hence, from this plot one expects three different phases 
each one corresponding to different density regions. 

Panels b) and c) of Fig. [1] are built from the PDFs 
shown in panel a). There, the dimensionless pressure 
Z = j3P/ p is plotted (dark circles) as well as the isother- 
mal compressibility Xi both as a function of the most 
frequent (p. The x values are obtained by the density 
fluctuations, i. e., by x = ^(< p^ > - < p >^)/ < p 
which should equal x — ^p/^{PP) according to the 
fluctuation-dissipation theorem. Finally, panel d) shows 
(dark circles) the order parameter Qe, as defined in 
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refs. [23I, HH, [26| , implemented for the geometric centers of 
the ehipsoids. The isotropic-nematic (hquid-hquid crys- 
tal) phase transition is observed as a plateau of Z{(p) 
and no discontinuity is detected. However, this transi- 
tion is clearly pointed out by the large jump of x{^)^ 
which peaks dit (p = (pi-n ^ 0.341. Additionally, Qei^) 
shows no sign of change at this density, indicating the 
expected absence of positional order. The second transi- 
tion occurs in the interval = 0.584 < cp < cps = 0.605 
and corresponds to Z ^ 27.0 (for TV ^ 00 we expect 
a shift of (ff^ Lpsi and Z to larger values [1^). It is 
characterized by a discontinuity of Z{(f)^ and a jump 
of x(^) and QQ^if). This Qq{^) jump is indicative of 
the appearance of positional order of the ellipsoids' geo- 
metric centers. It should be noted that the fluid-crystal 
transition interval for hard ellipsoids is well shifted to 
the right as compared to the hard spheres transition 
{(pf = 0.492 <ip <ips= 0.545) [13, [23, [Is^ . Something 
similar happens to the maximum hard ellipsoid packing 
fraction [29]. A fit of the form Z ~ {ipd — ^)~^ to the 
high pressure branch of Z{ip) leads to a divergence at 
ifd = 0.768. This value far exceeds the largest packing 
fraction of hard spheres, cp = 0.74, and is close to the 
maximum density of hard ellipsoids reported by Donev 
et. al., 0.77 [H. 

Panels b), c), and d) of Fig. [T] also include the re- 
sults obtained for the BP analytical approximation as 
light squares (these results are 3 times faster to obtain 
than those corresponding to the exact PW solution). As 
it can be seen, for 0.37 ^ (p ^ 0.57, the three plots, 
Z{(p)^ x{^)^ and Qei^) show a very good match be- 
tween BP and PW. However, both BP transitions are 
shifted to lower densities, which is a direct consequence of 
the systematical overestimation of the contact distance, 
c"BP ^ o-pw- Since BP clearly overestimates the T-shape 
configuration the equations of states quantitatively dis- 
agree for (p < 0.37, when T-shape configurations are rel- 
evant. This disagreement is clearly seen in panels b) and 
c) of Fig. [TJ Conversely, the T-shape configurations are 
practically absent in the nematic phase and thus the BP 
expression matches the exact PW solution. At very high 
pressures the systematic overestimation of a^p produces 
a less dense crystal branch, an underestimation of x{^)^ 
and a more structured crystal (see Qq{(p)). 

The structure of the different phases is analyzed in 
Fig. [2] where the radial distribution function and order 
parameter, p{r) =< l/2(3(ui • Uj)^ — 1) >, are shown to- 
gether with the corresponding snapshots for the highest 
(panel c)), the lowest a), and an intermediate pressure 
b). The corresponding phases are crystal, isotropic, and 
nematic, respectively. The isotropic phase is character- 
ized by a slowly increasing g{r) which peaks at r ^ a^. 
The p{r) shows the typical peak at r = (7\\ , which prac- 
tically vanishes at r ^ <Jx. Thus, there is not a long 
range alignment of the entities for cp < (pi-n- The corre- 
sponding snapshot clearly shows this fact (inset of panel 
a)). Conversely, the fluid structure for cpi-n < (p < (pf 



does show a long range alignment of the ellipsoids, as 
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FIG. 2. Radial distribution functions, g{r), (dark lines) and 
radial order parameters, p(r), (light lines) for a) the lowest, b) 
an intermediate, and c) the highest pressure. The insets are 
corresponding snapshots, d) Front view and, e), side view of 
the corresponding geometric centers of the inset at c) . Results 
are obtained for the exact apw- 

indicated by the large value of p{r) for all distances, r 
(panel b)). However, the ellipsoids' geometric centers 
still present a fluid-like structure, as pointed out by the 
g{r) function. Finally, for (p > (ps the g{r) grows sev- 
eral peaks pointing out a highly structured phase (panel 
c)). The first and larger one appears at r ^ 1.2, and 
corresponds to the touching and stacked ellipsoids shown 
in the inset. The second, at r ^ 2.4, is also due to the 
stacked arrangement and corresponds to two ellipsoids 
separated by a third one and sandwiched by them. The 
wide peak at r ^ 4.2 corresponds to side-to-side con- 
figurations of particles belonging to different stacks. It 
is located at r < ax since the particles of the adjacent 
stacks are partially sandwiched. The p{r) function also 
shows for this case the expected long range alignment. 
More details on the obtained crystal structure are seen 
in panels d) and e). Panel d) highlights the obtained 
hexagonal arrangement of the particles' geometric cen- 
ters (in correspondence to the snapshot inserted in panel 
c)). Additionally, panel e) shows the side view of the 
stacks' axes, which align defining planes. 

In summary, this work provides the equation of state 
of hard 1:5 aspect-ratio oblate ellipsoids by means of 
replica exchange simulations for a moderate system size 
{N = 100). This is done by considering the BP analyti- 
cal expression and the exact numerical solution for deter- 
mining overlaps. We observed a good overall agreement 
between BP and the exact numerical solution, though 
deviations appear at low concentrations, where T-shape 
configurations are frequently present in the system struc- 
ture. For both cases, an isotropic-nematic and a nematic- 
crystal transitions are captured, occurring for the exact 
case at (pi-n = 0.341 and Z = 8.2, and in the interval 
(pf - (Ps = 0.584 - 0.605 and Z = 27.0, respectively. 
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